Macroscopic consequences of calcium signaling in microdomains: 
A first passage time approach 
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Calcium (Ca) plays an important role in regulating various cellular processes. In a variety of 
cell types, Ca signaling occurs within microdomains where channels deliver localized pulses of Ca 
which activate a nearby collection of Ca-sensitive receptors. The small number of channels involved 
ensures that the signaling process is stochastic. The aggregate response of several thousand of 
these microdomains yields a whole-cell response which dictates the cell behavior. Here, we study 
analytically the statistical properties of a population of these microdomains in response to a trigger 
signal. We apply these results to understand the relationship between Ca influx and Ca release in 
cardiac cells. In this context, we use a first passage time approach to show analytically how Ca 
release in the whole cell depends on the single channel kinetics of Ca channels and the properties of 
microdomains. Using these results, we explain the underlying mechanism for the graded relationship 
between Ca influx and Ca release in cardiac cells. 



The use of Ca-sensitive receptors is ubiquitous in the 
design of signal transduction processes. A basic feature 
of this design is the close positioning of Ca-permeable 
channels to Ca-sensitive receptors in localized regions of 
the cell, often referred to as Ca microdomains [1]. Typi- 
cally, the signaling takes places when a "trigger" channel 
delivers Ca into the microdomain, and the rise of the lo- 
cal Ca concentration is sensed by Ca-sensitive receptor 
channels. In this way, a trigger input can induce a large 
response signal which then facilitates downstream cellu- 
lar processes. In cardiac cells this signaling architecture is 
used to mediate the interaction between membrane volt- 
age and cell contraction [2] , while in neurons it modulates 
synaptic function [3[. 

Within a microdomain the relationship between the 
trigger signal and receptor response is nonlinear, since 
(i) Ca receptors themselves control the flow of a high 
concentration of Ca sequestered inside the cell, and (ii) 
Ca receptors are clustered so that Ca released from one 
receptor can induce a neighboring receptor within the 
microdomain to open. In this way a small increase in 
trigger current can induce a large autocatalytic response 
from a cluster of receptors, leading to an "all or none" 
excitability. However, the spatial separation between mi- 
crodomains ensures that the Ca released in one domain 
does not necessarily induce release in a neighboring do- 
main. Thus, the cell's response to a signal such as a 
change in membrane potential will be dictated by the 
aggregate response of a population of independent mi- 
crodomains. This feature of Ca signaling, referred to in 
cardiac physiology as "local control" |4|, [5| , allows the cell 
to smoothly control its whole-cell response, despite the 
all-or-none nature of Ca signaling at the channel level. 

Detailed simulation studies have shown that the en- 
semble behavior of a population of microdomains can re- 
produce experiment ally- measured whole-cell currents [6, 




FIG. 1: A. Illustration of Ca- mediated signaling between LCC 
and RyR channels in cardiac cells. Ca is injected into the cell 
via the LCCs and triggers the opening of RyR channels in the 
immediate vicinity. B. Markov state models for LCC (left) 
and RyRs (right). The vertical dots indicate many (~ 50 — 
300) RyRs; horizontal dots represent possible deeper states. 

LZ|. However, these studies did not provide a concise an- 
alytic relationship between channel statistics and whole- 
cell behavior. Also, Ca receptor clusters have been stud- 
ied extensively but without directly address- 
ing the relationship between trigger and response. In 
this paper we use a first passage time approach to study 
the nonlinear stochastic interaction of Ca-sensitive re- 
ceptor channels and nearby "trigger" channels. We then 
compute the aggregate response of a population of mi- 
crodomains to show how the whole-cell response depends 
on local Ca signaling. We apply these results to explain 
experimentally known features of Ca signaling in cardiac 
cells. For the sake of concreteness we have developed 
our analysis within the context of Ca signaling in cardiac 
cells. However, the signaling architecture we consider is 
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found in many cellular systems [1] and should be appli- 
cable in a broader context. 

The basic architecture of Ca signaling in cardiac cells 
is shown schematically in Fig.[T]/L Here, the signaling oc- 
curs between voltage-sensitive L-type Ca channels (LCC) 
on the cell membrane and a cluster of Ca-sensitive Ryan- 
odine receptors (RyR) [2|. The RyRs gate the flow of Ca 
from an internal store, the sarcoplasmic reticulum (SR), 
which has a concentration roughly 10 4 times greater than 
in the cell. Within the microdomain, signaling occurs via 
the interaction of a few LCCs and a tightly knit cluster 
of 50 to 300 RyRs. The microdomain is shaped like a 
pill box of height ~ 10 nm and diameter ~ 100 nm; a 
typical cardiac cell has ~ 10 4 microdomains distributed 
throughout the cell. 

To model ion channels within the microdomain we 
use a Markov state approach as illustrated in Fig. [TJ>. 
The single channel properties of LCCs have been stud- 
ied extensively [ll|; here we model a single LCC using 
two closed states which can transition to a Ca-permeable 
open state . The membrane voltage (V) dependence 
of LCCs appears in the transition rates between closed 
states (ai(V), f3i(V)), while transition rates to and from 
the open state (a and (3) are voltage independent. The 
activation kinetics of RyRs are also well known [12| , and 
it is believed that the open probability is regulated by 
several Ca binding sites acting cooperatively on the re- 
ceptor. We use a minimal model which incorporates these 
essential features by making the open rate k+c 2 (where c 
denotes the local Ca concentration) and closing rate k- . 
The interaction between LCCs and RyRs is dominated 
by the activation kinetics, and so it is safe to neglect 
deeper states in the Markov scheme that describe slower 
inactivation processes. 

When a cardiac cell is stimulated the rise of the mem- 
brane voltage leads to a dramatic increase in the open 
probability of LCCs. The subsequent Ca entry into the 
cell induces Ca release via RyR receptors, leading to a 
global rise of Ca in the cell. High resolution optical imag- 
ing of Ca [13] in cardiac cells indeed reveals local eleva- 
tions in Ca, referred to as Ca "sparks". These sparks 
represent the release of Ca from the SR through RyR 
clusters which are activated by local LCCs. Ca sparks 
have relatively constant amplitude and duration [13], and 
so the release flux of Ca from the SR is largely dictated 
by the number of Ca sparks recruited in the cell [14(. If 
we denote ANs to be the number of new sparks recruited 
during a small time interval At, then 



AN S = N d P Cl aAtP s 



(1) 



where Nd is the number of microdomains in the cell; 
Pd^At is the probability that an LCC makes a tran- 
sition from C\ to O in the time interval At, with Pc 1 
being the probability of the closed state that directly 
transitions to the open state of an LCC; and Ps de- 
notes the probability that the LCC opening will trigger 



a spark. Thus, the key quantity relating the local chan- 
nel interactions to the whole-cell response is the rate of 
spark recruitment Rs = ANs/ At, which depends on the 
trigger-response interaction described by Ps. 

To compute Ps we first note that Ca ions inside the 
cell diffuse at ~ 150 — 300 /im 2 /s, and thus equilibrate 
over the scale of the microdomain much faster (~ 0.01ms) 
than the typical channel transition times (~ 1ms). Fast 
diffusion yields c w c Q + (r /v) (n • ir v r + k • ica), where 
c Q is the resting Ca concentration outside of the mi- 
crodomain of volume v; r is the time constant of diffusion 
out of the microdomain; n and k are the number of open 
RyR and LCC channels, respectively; and ir v r and ic a 
are the fluxes, in ions per unit time, through the respec- 
tive single open channels. We set ir v r = gc sri valid for 
small c, where g is the single RyR channel conductance 
and c sr is the Ca concentration in the SR. Since there are 
only a few LCC channels in the microdomain, each with 
small open probability [ll[, we wm simplify the system 
further and assume that there is only one LCC channel 
in each microdomain so that k = or 1. 

The state of the RyR cluster can be described using a 
master equation for the probability P(n, t) that n out of 
N total RyR channels in the cluster are open, given by 



dP(n, t) 
dt 



r + (n - l)P(n - 1, t) + r_ (n + l)P(n + 1, t) 



[r + (n) + r_ (n)] P(n, t) 



(2) 



with forward transition rate r + (n) = k+(N — n)c 2 (a 
third-order polynomial in n) and backward transition 
rate r_(n) = k-n. For N large we can use well-known 
approximation methods [15[ to reduce this one-step pro- 
cess to its corresponding Fokker-Planck equation (FPE), 
which we can write, in terms of the fraction x = n/N of 
open channels, as 



dp(x,t) 
dt 



d_ 

dx 



[f(x)p(x,i)] 



1 d 2 

2N dx 2 



[h(x)p(x,t)],(3) 



with drift comcient f(x) = (r + (7Vx) - r_(7Vx))/7V and 
diffusion coefficient h(x)/N = (r+(Nx) + r-(Nx))/N 2 . 
The diffusion term scales as 1/N implying the increas- 
ingly deterministic nature with larger cluster size [16]. 

The equivalent Langevin description corresponding to 
this FPE implies time evolution of the mean open frac- 
tion (x) dictated by the equation = f((x)). There- 
fore f(x) = determines the stationary points of x. 
The deterministic component of the dynamics can be 
conveniently visualized in terms of a potential U(x) = 
-j*f(x')dx' (Fig. El). In the case when the Ca chan- 
nel is closed (k = 0), the potential landscape reveals 
a bistable system with three fixed points. The stable 
point x a corresponds to an inactive cluster with very 
few channels open. A second stable point x c denotes 
a nearly fully open cluster. Finally, an unstable point 
at position x^ indicates the peak of a potential barrier 
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FIG. 2: Main: Effective potential landscape U(x) correspond- 
ing to closed (k = 0) or open (k — 1) LCC. Inset: Barrier 
near origin disappears when ic a is large enough. Parameters 
here are ic a = 0.06 pA, g — 0.520 fim 3 /is -1 and c Q = 5 \iM. 



separating the closed and open cluster states. Since x 
is small during initial spark activation, to second or- 
der we have f(x) ~ a + \ix + k+q 2 x 2 where a = k+s 2 
and \i = — k- + 2k+qs, with g = N(r/v)gc sr and 
5 = c + (r/v)kica- The equilibrium states can then 
be found as x a « (k + /k-)c 2 and £5 « (k_/k+)q 2 . 

The interaction between the LCC and the RyR clus- 
ter is governed by the degree to which an LCC opening 
will induce a transition from x a to x b . When the LCC 
channel is open the barrier is reduced (or eliminated) as 
shown in Fig. [2j owing to the elevated Ca concentration 
in the microdomain due to Ca influx via the open LCC. 
If the LCC remains open until x crosses x bl then a spark 
will almost certainly occur as x proceeds to x c . On the 
other hand if the channel closes before such crossing has 
occurred, then as the barrier returns the cluster state will 
most likely return to x a . Thus, given a passage time t 
from x a to x bl the LCC open time, which follows an ex- 
ponential distribution fo(t') = /3exp(— /3t f ), must exceed 
t; this yields 



poo 

P s ((3,x a ,x b ) = / dt e~ (3t P(x a ,x b ;t), 
Jo 



(4) 



where P(x a ,x b ;t) is the first passage time density 
(FPTD) from x a to x b for the system given by Eq. ([3j). 

The statistics of this FPTD have been worked out in 
detail in a classic paper by Darling and Siegert [17[ , where 
they show that the Laplace transform of the FPTD can 
be found by solving the adjoint FPE. Linearizing the 
drift and diffusion coefficients so that f(x) = a + \ix 
and h(x) = a + 7#, where 7 = fc_ + 2k+qs, the ad- 
joint FPE can be written as a confluent hypergeomet- 
ric differential equation, which has the linearly indepen- 
dent solutions u\(x) = M(mi : rri2 : y(x)) and U2(x) — 
y(x) 1 1712 M(l-\-mi — 77i2, 2 — 772,2, y( x )), where M denotes 
Kummer's function [18[ with arguments mi = — /3//i, 
ra 2 = 2Na(j — /-O/7 2 , and y(x) = —2N/i(a + jx)/j 2 . 
Imposing reflecting boundary conditions at x = 0, since 
the master equation (Eq. (|2J)) does not admit negative 
valued states, the probability of sparking can be writ- 
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FIG. 3: Analytical (curves) and numerically simulated (sym- 
bols) sparking probability Ps as a function of ica , conditional 
upon an LCC activation, for conductances g — 0.910 (■), 
0.628 (□), 0.491 (A), and 0.416 (A) /im 3 /^" 1 - 



ten [17[ as 



P s (P,x a ,x b ) 



Ui(Xa)/u[(0) ~ U 2 {x a )/u' 2 {0) 

u^x^/u^O) - u 2 (x b )/u' 2 (0) ' 



(5) 



where 1^(0) and u' 2 (0) are the derivatives of the inde- 
pendent solutions evaluated at the origin. Eq. ([5j) gives a 
full description of the sparking probability given a trigger 
current ic a with duration from the distribution fo- 

The analytic formula for the spark rate (Eq. (|5])) can 
be simplified further in the limit of large or small LCC 
current (ic a ) . Using the asymptotic form of the Kummer 
function [18[, for large %c a we have 



Ps 



cr + jx b J 



(6) 



In this limit, the response is governed by the mean open 
time of the LCC, and does not strongly depend on the 
stochastic properties of the RyR cluster. For small ic a a 
barrier must be surmounted, and stochastic fluctuations 
due to the N RyR channels dictate the probability for 
sparking. In this regime, to leading order we have 



P s - exp [—N(x b - x a ) 
V 7 



(7) 



showing the strong dependence on RyR cluster size N. 

In order to check the validity of this analytical result 
we simulate independent Monte Carlo trajectories of the 
master equation to estimate the probability of sparking 
(Ps) given an LCC opening at time t = 0, in an en- 
semble of 100, 000 microdomains whose clusters are ini- 
tially closed (x = x a ). Within each unit a single LCC is 
opened for a duration t' chosen from fo- We then esti- 
mate Ps as the fraction of microdomains whose clusters 
reach n = x b N within that time. In Fig.[3]the predictions 
of Eq. ([5]) are compared with the numerical simulations. 
As shown, the agreement is quite good. 

In experiments, the relationship between Ca entry and 
Ca release can be assessed by depolarizing the membrane 



4 




FIG. 4: Main: Analytical (left curve) and numerically sim- 
ulated (■) spark recruitment rate Rs as a function of Vmax 
(g — 0.910 /im 3 /xs _1 ); whole-cell I ca (right curve). Inset: Ex- 
perimental values from [5[ for peak SR release current (•) and 
peak LCC current (o). (All plots normalized to peak height) 



level. The main result is an analytic description of the 
relationship between single channel kinetics and the ag- 
gregate whole-cell response. In the context of the cardiac 
cell, we have reproduced important voltage-dependent 
relationships which are known experimentally, but un- 
til now have not been explained in detail. This work 
should pave the way to a more detailed understanding of 
Ca signaling in a wide range of biological processes. 

This work is supported by the NSF and NIH/NHLBI 
P01 HL078931; R.R. thanks Tom Chou and Vladimir 
Minin for valuable discussions. Y.S. thanks A. Karma 
for valuable discussions, and the KITP Santa Barbara, 
where part of this work was completed. 



to various test voltages V max . The peak total currents for 
LCC Ca entry and SR Ca release are then measured fol- 
lowing the depolarization. In Fig. 4 (inset) we reproduce 
the experimental data of Wier et al. [5j], showing that 
both Ca release and Ca entry follow a bell-shaped de- 
pendence on Vmax- This relationship is well established, 
and referred to as graded release. Since the spark re- 
cruitment rate is directly related to the amount of Ca 
released from the SR, we check this relationship by com- 
puting the spark rate Rs as a function of Vmax- Figure 
[4] shows the peak spark rates due to a depolarization to 
Vmax , using our analytically calculated Ps from Eq. ([5j) in 
Eq. (pQ), as well as by simulating Nd = 100,000 indepen- 
dent microdomains with Markov-governed LCCs initially 
in state C2 • Numerically, the spark rates were estimated 
by counting the number of sparks recruited in 1ms in- 
tervals. Also plotted is the peak whole-cell Ca current 
estimated by Ic a = NdP ica, where P Q is the steady- 
state open probability of the LCC. As shown, consistent 
with graded release, the spark rate dependence on V max 
follows a bell curve mirroring the Ca entry Ica, an d the 
analytical predictions agree semi-quantitatively with the 
numerical simulations. 

The graded relationship between whole-cell Ca release 
and Ca entry is reflected in the spark recruitment rate 
(Eq. JT])). For large negative Vmax ica is large and every 
LCC opening will trigger a spark. In this regime the 
spark rate is dominated by the voltage dependence of 
Pci, which roughly follows P G , since a and (3 are fast 
voltage-independent rates. As V ma x is increased further, 
ica decreases and Ps is reduced as the probability of 
triggering a spark is less. Similarly, Ic a also decreases 
due to the drop in the single channel current ic a - Hence, 
Eq. (pQ) reproduces experimentally measured whole-cell 
properties in terms of the kinetics of ion channels in Ca 
microdomains. A more detailed analysis of the analytic 
prediction will be presented in a future publication. 

In this letter we have analyzed the stochastic proper- 
ties of Ca signaling at both the whole-cell and ion channel 
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